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NOTATIONS 


E 

total  energy  per  unit  of  mass  (J/kg) 

Tf 

propellant  flame  temperature  (K) 

F 

convective  flux  vector  in  the  x  direction 

u 

velocity  component  in  the  x-longitudinal 

G 

convective  flux  vector  in  the  y  direction 

direction  (m/s) 

H 

convective  flux  vector  in  the  z  direction 

Uy 

axial  mean  velocity  (m/s) 

S 

source  term  vector 

V 

velocity  component  in  the  v- lateral 
direction  (m/s) 

k 

turbulent  kinetic  energy 

Vinj 

wall  injection  velocity  (m/s) 

< 

mass  flow  rate  per  injecting  surface  unit 
(kg/m2) 

W 

velocity  component  in  the  iv-lateral 
direction  (m/s) 

P 

pressure  (Pa) 

x,y,z 

co-ordinate  system  (m) 

P’ 

Q 

fluctuating  pressure  (Pa) 

conservative  variables  vector  of  Navier- 

8 

turbident  dissipation  rate 

Stokes  equations 

K 

thermal  conductivity  (W/mK) 

m 

total  mass  flow  rate  (kg/s) 

P 

dynamic  viscosity  (kg/ms) 

t'b 

burning  rate  (mm/s) 

P 

density  (kg/m3) 

T 

temperature  of  the  flow  (K) 

Ps 

propellant  density  (kg/m3) 

Dimensionless  Parameter 

Rec  axial  Reynolds  number:  p  uv  DJp,  where  Dc  is  a  characteristic  diameter  of  the  SRM 
Res  wall  injection  Reynolds  number:  p  vw  Dc/p 

y  isentropic  exponent,  equal  to  the  specific  heat  ratio  for  a  perfect  gas 

INTRODUCTION  [1] 

Internal  ballistics  in  a  SRM  can  be  solved  with  various  ways  and  for  various  objectives  [2].  The  motor 
design  engineer  wants  to  predict  or  understand  the  burning  characteristics  and  the  global  performances  of 
the  motor,  seek  the  efficiency  of  thermal  insulation  and  nozzle  design,  check  the  reliability  of  ignition  and 
motor  design  for  the  life  cycle  of  the  motor,  take  into  account  variability  induced  by  manufacturing 
processes,  etc. 


Paper  presented  at  the  RTO/VK1  Special  Course  on  “ Internal  Aerodynamics  in  Solid  Rocket  Propulsion  ", 
held  in  Rhode-Saint-Genese,  Belgium,  27-31  May  2002,  and  published  in  RTO-EN-023. 
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ORGANIZATION 


Different  ways  can  be  followed  to  reach  the  same  goal.  The  reasons  why  a  project  manager  uses  numerical 
simulations  are  multiple.  In  the  predictive  mode,  the  aim  is  to  design  the  system  which  fulfils  the 
specifications  at  the  lowest  development  cost  (by  minimizing  the  number  of  prototypes  and  tests).  In  the 
explanatory  mode,  the  desire  is  to  explain  an  observed  and  unknown  phenomenon  that  happened  during 
the  development  phase  of  a  program. 

For  solid  rocket  motors,  the  objectives  of  internal  flow  computations  are  to  predict  the  global 
performance  and  reliability  of  the  rocket,  avoid,  minimize  or  master  undesired  behaviors  (like  thrust 
oscillations,  erosive  burning),  take  into  account  flow/grain  and  casing  interactions  (mechanical  and 
thermal  loads),  etc. 

The  internal  aerodynamics  inside  a  solid  rocket  motor  can  be  modeled  with  increasing  degrees  of 
complexity  from  the  simplest  global  equations  to  a  full  3D  numerical  simulations.  An  AGARD  Lecture 
Series  has  already  been  organized  on  the  Design  Methods  in  Solid  Rocket  Motors  [3].  The  objective  of 
this  special  course  is  not  to  duplicate  the  materials  developed  in  AGARD-LS-150.  We  will  focus  on 
multidimensional  simulations  of  internal  flows. 

Example  of  objectives  and  constraints: 

•  assessment  of  stability 

•  prediction  of  performances 

•  prediction  of  reliability 

•  prediction  of  variability,  for  instance  thrust  imbalance  (important  when  using  simultaneously 
more  than  one  identical  motor) 

Conception: 

•  choose  the  right  propellant  for  the  application,  and  design  the  initial  geometry  of  the  solid 
propellant  charge  that  will  deliver,  with  surface  regression,  the  required  time  history  of  the  gas 
flow 

•  assess  reliability  of  the  designed  motor  (mechanical  loads,  ignition,  casting  process  and  raw 
materials  variability,  thrust  oscillations,  . . .) 

Challenges  of  numerical  modeling: 

•  two-phase  reacting  flow  (aluminized  propellants) 

•  gas  in  a  wide  range  of  temperatures 

•  multi-species  and  turbulence 

•  two  or  three  dimensional  geometry 

•  moving  boundaries 

•  fluid/structure  coupling,  with  heterogeneous  surface  combustion 

•  steady  and  unsteady  compressible  flows,  with  all  range  of  Mach  number 

The  impact  of  numerical  modeling  in  designing  motors  of  new  generation  had  a  recent  increase  due  to 
mainly  three  causes: 

•  large  progress  of  computer  power  (hardware)  and  in  computational  fluid  dynamics  (software) 
in  the  last  decades 

•  objective  of  cost  reduction  (mainly  in  the  development  phase  of  a  new  motor) 

•  new  objectives  of  performance  and  reliability 
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An  improved  prediction  reduces  the  number  of  qualifying  fires,  and  for  large  and  expensive  motors,  this  is 
of  prime  importance. 

This  lecture  series  will  be  divided  in  two  special  courses.  The  first  one  presents  the  general  models  for 
solving  internal  steady  state  aerodynamic  in  solid  rocket  motors,  the  second  one  focus  on  pressure  and 
thrust  oscillations  modeling. 


GENERAL  EQUATION  FOR  AERODYNAMICS 

At  a  starting  point,  the  general  equations  describing  fluid  flows  inside  a  solid  rocket  motor  begins  with  the 
Navier-Stokes  equation.  We  remind  them  as  an  introduction. 

The  general  form  of  the  conservation  equations  for  a  three  dimensional  viscous  flow  can  be  written  in  the 
following  form: 

6Q  8F  dG  dH  0 

—  +  —  +  —  + - =  S 

dt  dx  By  dz 


where: 


Q  = 


p 

pu 

pv 

pw 


pE 


pii 

pir  +P-crxx 


F  = 


puv-aXJ 


pUW  —  <7  xz 

dT 

(pE  +  P)u  -  uaxx  -  V(Txv  -  wcrxz  ~k  — 
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pv 
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H  = 


pw 

pUW-C Txz 

pvw  -  a  yz 
pw2  +P-cjz: 


{pE  +  P)w-uctx:  - vay: 


-  WCT 


zz 


dT 

—  K - 

dz . 


and  S  is  the  source  term  vector. 

These  equations  must  be  solved  with  appropriate  solvers  and  on  appropriate  grids. 

Over  the  last  twenty  years,  different  grid  technologies  have  been  developed.  Based  on  experience  on  finite 
difference  schemes  and  for  the  sake  of  simplicity  in  software  management,  structured  grids  have  been 
widely  used.  A  structured  grid  is  a  grid  where  each  component  can  be  identified  by  two  indices  (i,j)  in  2D 
and  three  indices  (i,j,k)  in  3D.  Each  index  refers  to  a  space  dimension.  Hence,  the  grid  points  are  ordered, 
giving  the  name  of  “structured  grid”.  However,  this  technology  was  shown  very  precise  and  efficient  for  a 
simple  geometry  but  limited  for  a  complex  geometry.  Techniques  of  multiple  structured  domains 
overlapping  have  been  developed,  but  with  special  difficulties  for  ensuring  global  conservation  among  the 
different  domains  and  a  good  efficiency  on  distributed  memory  parallel  computers. 

For  an  unstructured  flow  solver,  the  computational  domain  is  tessellated  using  a  grid  composed  of 
simplices,  which  are  quadrilaterals  or  triangles  in  two  dimensions  and  generally  tetrahedras,  pyramids, 
pentagons,  prisms  and  hexahedras  in  three  dimensions.  Unstructured  grids  provide  flexibility  for 
tessellating  about  complex  geometry  and  for  adapting  to  flow  features,  such  as  shocks  and  boundary 
layers. 

On  a  given  grid,  one  has  the  option  of  locating  the  variables  at  the  cell  centers  or  at  the  vertices  of  the  grid, 
giving  rise  to  cell-centered  and  cell-vertex  schemes.  Alternatively,  it  is  possible  to  deal  strictly  with 
averages  defined  over  volumes.  This  approach  has  certain  advantages  for  higher  order  schemes.  In  the 
case  of  finite  volume  schemes,  the  governing  equations  are  discretized.  This  allows  discontinuities  to  be 
captured  as  part  of  the  solution. 


SOLID  PROPULSION  MODELS 

General  modem  CFD  codes  for  computing  flows  in  solid  rocket  motors  have  the  following  features: 

•  solves  the  2D  axisymmetrical,  plane  and  3D  Navier-Stokes  equations  for  laminar  or  turbulent 
flows 

•  uses  unstructured  meshes  for  complex  geometry  treatment 

•  has  the  possibility  for  treating  the  chemical  reactions  of  multi-species  and  the  coupling  between  a 
gas  phase  and  a  condensed  phase,  inert  or  not,  with  specific  models 

•  has  moving  mesh  facilities 

•  incorporates  specific  solid  propulsion  models  for  the  burning  rate,  from  simple  laws  (regression 
rate)  to  several  coupling  (ignition,  erosive  burning,  unsteady  combustion)  as  well  as  solid 
propellant  grain  coupling  (mechanical,  surface-bumback). 


9-4 


RTO-EN-023 


Numerical  Modeling  of  Internal  Flow  Aerodynamics 
Part  1:  Steady  State  Computations 


COMBUSTION  AND  EROSIVE  BURNING 

Erosive  burning  is  a  phenomenon  commonly  experienced  in  a  solid  propellant  rocket  motor,  represented 
by  an  increase  of  the  local  propellant  burning  rate  due  to  high  velocity  combustion  gas  flow  across  the 
burning  surface.  Most  propellants  have  a  minimum  cross-flow  velocity  below  which  erosive  burning  is  not 
observed,  referred  to  as  the  “threshold  velocity”. 

The  erosive  burning  mechanism  is  believed  to  be  due  to  the: 

•  increase  in  gas-to-solid  heat  feedback  caused  by  the  increase  in  transport  coefficients 

*  turbulence-enhanced  mixing  and  chemical  reaction  of  the  oxidizer  and  fuel  rich  gases  pyrolized 
for  composite  propellants 

Steady  combustion  is  a  complex  mechanism  including  chemical  and  physical  effects  (nature  and  details 
of  energetic  materials  and  additives,  particle  size  distribution,  operating  conditions:  pressure,  initial 
temperature,  radiation,  ...).  For  erosive  burning,  the  cross-flow  velocity  (parallel  to  the  solid  propellant 
burning  surface)  constitutes  an  additional  operating  condition  of  extreme  importance. 

Several  theoretical  approaches  have  been  reported,  which  can  be  grouped  in  five  categories,  following 
Kuo  and  coworkers  [4]: 

1)  phenomenological  heat  transfer  theories 

2)  modification  of  the  propellant  combustion  mechanism 

3)  integral  boundary  layer  analysis 

4)  chemically  reacting  turbulent  boundary  layer  analysis 

5)  others 

We  focus  on  models  1  and  4.  Models  in  the  category  1  will  have  an  interest  for  engineering  design 
problems  while  models  in  the  category  4,  relying  on  more  fundamental  viewpoint,  are  thought  to  be  more 
precise  and  thus  more  appropriate  for  being  incorporated  in  complex  CFD  codes. 


Phenomenological  Heat  Transfer  Theories 

Most  of  these  models  are  based  on  or  derived  from  Lenoir  &  Robillard  [5]  approach,  giving  the  general 
expression  for  the  burning  rate  in  the  form: 


=  aPn  + 


aG8 

- e  G 

m,) 


where  a,  /3  and  g  are  three  constants,  G  is  the  mass  flux  through  the  port,  and  Dh  is  the  hydraulic  diameter. 
The  expression  f(Dh)  can  include  scale  effects.  In  the  original  form,  f(Dh)=Dh0'2,  j3=53  is  found  to  be 
independent  of  the  propellant  type,  and  g=0.8  (based  on  Chilton-Cobum  correlation  for  evaluating  the 
convective  heat  transfer  coefficient ). 


This  approach  is  well  suited  for  ID  analysis. 


Chemically  Reacting  Turbulent  Boundary  Layer  Analysis 

These  models  are  well  suited  to  aerothermochemical  analysis  of  erosive  burning  of  composite  propellants, 
and  many  authors  have  contributed  (King  [6],  Beddini,  Kuo,  and  people  from  ONERA).  We  will  focus  on 
ONERA  approach  [7,8]  used  at  SNPE. 
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The  turbulent  flow  is  solved  in  a  high  Reynolds  approach  and  the  propellant  surface  is  treated  as  a  wall 
zone.  In  this  region,  Couette  equations  are  solved,  from  the  propellant  surface  to  the  first  integration  point 
of  the  flow  in  the  port.  The  temperature  gradient  in  the  flame  is  then  computed,  leading  to  the  heat  flux  to 
the  surface.  A  flame  height  criterion  is  used,  assuming  that  the  combustion  between  oxidizing  and  fuel 
gases  is  limited  by  the  diffusion  (valid  for  medium  and  large  AP  sizes).  In  the  erosive  regime,  the  solution 
of  the  Couette  flow  is  coupled  with  the  flame  height  criterion,  including  the  turbulent  contribution. 
This  coupled  system  is  solved  in  an  iterative  way,  until  the  velocity  and  temperature  profiles  match  the 
values  of  the  first  integration  point  in  the  flow.  At  convergence,  the  erosive  burning  rate  is  immediately 
obtained. 

Since  it  is  driven  by  viscous  effects,  the  erosive  burning  in  a  SRM  will  be  sensitive  to  the  scale  of  the 
motor. 


Temperature  T  (K) 


0  si/i 
0  si/io 

0  SI/IOO 

U  Calculated 

T  -  profiles 


a)  Scale  effect  on  the  flame  zone 


Figure  1:  Flame  Zone  Computed  by  the  ONERA  Model  at  Different  Scales  (1, 1/10  and  1/100). 


Figure  2:  Erosive  Behavior  as  a  Function  of  the  Motor  Scale  (ONERA  Model). 
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■ — ■ 

Slab  motors 

Marklund-Lake 

CT7-1 

Axisym  nozzleless  motor 

mmm 

Generator/sample  set-up 

Erosive  threshold  from  model 

c)  Erosive  mass  flow  rate  threshold  as  from  the  model 
and  various  small  scale  experiments 


Figure  3:  Erosive  Threshold  in  Various  Configurations  (ONERA  Model). 


TURBULENCE 

Turbulence  modeling  in  the  port  becomes  necessary  for  the  evaluation  of  heat  transfer  and/or  diffusion 
related  phenomena  (erosive  burning,  heat  flux  over  thermal  inhibitors  and  material  decomposition,  . . .). 

Cold  flow  experiments  with  wall  injection  have  shown  that  this  kind  of  flow  have  a  delayed  turbulent 
transition,  as  illustrated  on  Figure  4.  The  important  parameter  is  the  injection  Reynolds  number 
(defined  with  injection  velocity  at  the  blowing  surface  and  port  radius).  For  injection  Reynolds  number 
above  50,  the  transition  from  a  laminar  to  a  turbulent  flow  is  delayed. 


Figure  4:  Laminar  /  Turbulent  Transition  as  a  Function  of  the  Injection  Reynolds  Number. 

The  difficulty  in  simulating  turbulence  in  a  SRM  comes  from  the  fact  that  the  transition  is  always  inside 
the  port,  since  the  velocity  at  the  head-end  is  equal  to  zero.  So,  turbulence  modeling  must  compute 
correctly  the  transition.  A  schematic  view  of  laminar-turbulence  interaction  is  given  in  Figure  5. 
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Figure  5:  Schematic  View  of  Laminar-Turbulence  Transition  in  a  SRM. 


Usually,  classical  turbulence  models  are  used  for  internal  flows  in  SRM.  As  in  the  usual  treatment  of 
turbulence,  the  velocity  field  u  and  pressure  P  are  decomposed  into  mean  u,  P  and  fluctuating  u  \  P’  parts 
with  Favre’s  average  for  compressible  flows. 

Two-equation  models,  like  the  k-s  model,  solves  transport  equations  for  the  turbulent  kinetic  energy  k  and 
the  dissipation  rate  s.  Source  terms  (S)  are  accounted  for  modeling  turbulence  creation  and  dissipation. 
The  eddy  viscosity  is  expressed  as  a  function  of  k2/s  and  the  Boussinesq  hypothesis  is  used  for  computing 
the  turbulent  stresses. 

Classical  k-s  turbulence  models  are  isotropic,  and  this  hypothesis  is  very  restrictive  for  flows  in  SRM. 
Anisotropic  turbulence  models,  like  Algebraic  Stress  Models  (ASM)  are  an  efficient  way  for  improving 
turbulence  modeling  without  a  dramatic  CPU  increase. 

Flow  turbulence  modeling  must  always  be  associated  with  a  consistent  injection  modeling.  As  the 
difficulty  is  in  predicting  turbulence  transition,  the  flow  features  must  be  well  described  at  injection. 
This  can  be  done  by  solving  the  boundary  zone  with  Couette  or  Prantl  equations,  as  for  the  erosive 
burning  modeling,  or  with  approximated  laws  [9], 


Figure  6:  Computation  of  Turbulence  Level  in  a  Finocyl  SRM  (SNPE). 
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TWO-PHASE  FLOW  EFFECTS 

Many  solid  rocket  motors  use  aluminized  propellants  in  order  to  improve  their  performances. 
The  aluminum  combustion  produce  a  condensed  phase  (aluminum  oxide),  and  therefore  a  two-phase  flow 
in  the  rocket  chamber. 


An  important  consequence  of  the  presence  of  this  liquid  phase  is  two-phase  losses,  and  in  the  case  of  large 
segmented  motors  with  a  submerged  nozzle,  slag  accumulation,  which  may  have  several  consequences  on 
specific  impulse,  thermal  insulation  behavior  and  thrust  vectoring. 

Continuous  efforts  have  been  done,  for  many  years,  in  the  USA,  in  Russia  and  in  France,  in  order  to 
develop  numerical  models  able  to  predict  accurately  losses  and  the  slag  formation. 

More  details  can  be  fond  in  Salita  special  course  [10]. 


An  eulerian  or  a  lagrangian  description  can  be  used  for  the  condensed  phase. 


In  the  eulerian  description,  the  condensed  phase  is  assumed  to  be  a  continuous  medium,  and  the 
conservation  equation  are  derived  from  integrals  of  conservative  quantities  (mass,  momentum  and  energy) 
on  control  volumes.  The  following  vectors  are  added  to  the  gas  phase  vector  in  the  previous  form  of  the 
conservation  equation  [11]: 


QP  = 


pp 

Ppup 

Ppvp 


Ppwp 

Ppep 


PpUp 

1 

^5 

h3 

1 

Ppwp 

PpUp 1 

Ppupvp 

Ppupwp 

Ppupvp 

GP  = 

Ppvp2 

HP  = 

p  V  w 

r  p  p  p 

Ppupwp 

Ppvpwp 

PpWp 

(ppep\p_ 

(ppEp\p_ 

1 

Ja 

fcT 

3 

_ 1 

and  the  source  term  S  account  for  gas-condensed  phase  interactions.  In  these  expressions,  subscript  p  refers 
to  the  particulate  phase.  The  volume  fraction  of  the  dispersed  phase  is  noted  ap  and  is  supposed  to  be 
small  (for  being  neglected  in  the  gas  phase  equations).  Apparent  condensed  phase  density,  pp,  is  equal  to 
pp=appR  where  p^  is  the  condensed  phase  true  density. 

The  eulerian  method  is  well  suited  for  particles  with  a  fixed  diameter.  Generally,  aluminum  oxide  particle 
show  a  bi  or  trimodal  distribution  after  aluminum  combustion  (resulting  from  smoke,  nominal  aluminum 
residues  and  agglomerated  aluminum  residues).  According  to  its  small  size  and  relaxation  time,  the  smoke 
can  be  treated  with  the  equivalent  gas.  Generally,  one  or  two  classes  (diameter)  or  particles  are  used  in  the 
computation.  Even  if  it  is  possible  in  the  eulerian  form,  for  coupling  the  particle  velocity  field  with 
turbulence  or  taking  into  account  complex  phenomena  such  as  coalescence  or  break-up,  the  lagrangian 
method  is  more  appropriate. 
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In  the  lagrangian  method,  group  of  particles  are  emitted  from  the  propellant  surface  and  are  explicitly 
followed  in  the  flow.  It  allows  more  complex  physics  to  be  taken  into  account  since  it  approaches  the 
“discrete  form”  of  the  condensed  phase,  but  the  lagrangian  method  is  CPU  time  consuming  for  a  correct 
treatment  of  the  flow  (approaching  the  apparent  particle  density  in  the  flow). 


tcdon  1  Melon?  Melon  3  Mdon4 


Figure  7:  Lagrangian  Computation  of  the  Two-Phase  Flow  in  Ariane  5  SRM  [13]. 


Two-Phase  Losses 

These  losses  are  created  by  the  non  equilibrium  between  gas  and  condensed  phase  from  the  chamber  to  the 
throat.  In  the  converging  part  of  the  nozzle,  the  flow  is  accelerated  and  the  condensed  phase  is  accelerated 
by  the  flow.  An  important  parameter  will  be  the  ratio  of  condensed  phase  relaxation  time  over  the  transit 
time  in  the  nozzle.  If  we  note  R  this  ratio: 

•  for  small  R,  the  particles  are  in  equilibrium  with  the  gas,  they  are  accelerated  and  give  their 
thermal  energy  (temperature)  to  the  gas,  contributing  to  the  impulse; 

•  for  large  R,  particles  will  be  in  non-equilibrium  with  the  gas,  and  an  impulse  deficit  will  be 
created. 

Additional  phenomena  to  take  into  account  are  particles  break-up  and  phase  changes. 

Investigation  of  the  Slag  Formation  [12] 

Slag  is  generated,  either  when  aluminum  oxide  droplets  impinge  on  some  portion  of  the  back-face  of  the 
nozzle,  or  when  they  are  captured  in  the  recirculation  zone  behind  the  submerged  nozzle. 
This  phenomenon  generally  produces  a  slag  pool  in  the  aft  end  of  the  motor  and  an  aluminum  oxide  liquid 
film  on  some  part  of  the  nozzle  wall. 


Figure  8:  Schematic  View  of  Main  Phenomena  Leading  to  Slag  Formation. 
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Different  strategies  can  be  applied  to  the  slag  formation  investigation:  a  fully  coupled  numerical 

computation,  or  a  two  step,  uncoupled,  computation. 

The  principle  of  a  two-step  uncoupled  calculation  can  be  summarized  as  follows  [13]: 

1)  The  steady  state  of  the  gas  flowfield  must  be  calculated  with  an  appropriate  turbulence  model; 
the  presence  of  the  liquid  phase  must  be  taken  into  account  by  assuming  that  the  two  phases  are  at 
equilibrium  everywhere  in  the  booster  (same  temperature,  same  velocity  and  constant  density 
ratio  deduced  from  the  propellant  composition),  using  the  equivalent  gas  (see  chapter  on 
thermochemistry). 

2)  The  condensed  phase  consists  only  of  aluminum  oxide  droplets,  ejected  directly  from  the 
propellant  surface;  when  the  combustion  zone  is  small  compared  to  motor  geometry, 
no  combustion  model  need  to  be  used;  the  droplet  size  distributions,  used  in  the  simulations,  must 
follow  the  experimental  distribution,  measured  with  a  quench  bomb. 

3)  If  the  motion  of  the  droplets  is  simulated  by  a  Lagrangian  method,  all  the  trajectories  are 
calculated  on  the  steady  gas  flowfield  and  a  stochastic  model  can  be  used  to  take  into  account  the 
influence  of  the  turbulence  field  on  the  droplet  dispersion. 

4)  The  total  slag  rate  can  be  calculated  by  summing  the  weights  of  all  the  droplets  which  impinged 
the  nozzle  back-face  in  the  stagnation  zone  or  which  are  trapped  in  the  recirculation  zone  in  the  aft 
end  of  the  motor. 

If  the  motor  is  unstable,  a  coupling  may  exists  between  vortex-shedding  and  particles  behavior  in  the  flow. 

In  that  particular  case,  an  unsteady  computation  of  this  coupling  is  necessary  [14]. 


Figure  9:  Computation  of  the  Coupling  between  Vortex-Shedding  and  Two-Phase 
Flow  in  an  Unstable  Segmented  Solid  Rocket  Motor  (Particles  Volume  Fraction). 


THERMOCHEMISTRY 

Since  the  flame  zone  in  SRM  is  very  small  compared  to  motor  length,  we  can  consider  that  the  injected 
gas  from  the  burning  surface  are  the  final  combustion  products. 

Generally,  their  thermochemical  properties  are  computed  from  equilibrium  computer  codes,  most  of  them 
based  on  the  original  Gordon  and  Me  Bride  CEC71  code  [15].  Since  classical  computations  are  done  with 
constant  thermodynamic  properties,  the  validity  of  assumptions  made  for  computing  them  must  be 
checked.  Classically,  two  thermodynamic  parameters  are  taken:  the  frozen  heat  capacity  and  y  in  the 
motor.  As  an  illustration,  the  evolution  of  the  frozen  heat  capacity  of  a  composite  solid  propellant  is  given 
as  a  function  of  the  Mach  number  in  a  SRM  (from  the  chamber  to  the  throat)  in  Figure  10. 
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Figure  10:  Example  of  the  Evolution  of  the  Heat  Capacity  with 
the  Mach  Number  in  the  Convergent  Part  of  a  Nozzle. 


Another  approximation  can  be  made  by  writing  the  combination  of  Saint-Venant  and  enthalpy 
conservation  equation: 

Tt  =7Y1  +  ^— lM2) 
f  2 

and  using  this  equation  at  the  throat  (M=l)  for  computing  an  approximated  y  with  the  value  of  the 
temperature  at  the  throat  computed  by  the  thermodynamical  code. 

This  equation  also  shows  that  temperature  and  velocity  or  pressure  and  density  are  correlated  in  a  SRM, 
and  that  temperature  and  pressure  are  independant. 

The  other  thermodynamic  parameter,  for  instance  the  specific  heat  capacity,  can  be  adjusted  for  giving  the 
correct  characteristic  velocity  c*,  hence  the  correct  pressure. 


eps  (A/At) 


Figure  11:  Comparison  of  Theoretical  and  Approximated  Pressure  in  a  Nozzle. 
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Figure  12:  Comparison  of  Theoretical  and  Approximated  Pressure  in  a  Nozzle. 


For  a  two-phase  flow  approximation,  the  following  relations  must  be  written  for  the  heat  capacity  and 
specific  heat  ratio: 


=  (1  +  Cm)c-CmcB 


r 


g  _ 


1 


i-a+cj 


7~  1 
7 


where  the  subscripts  s  and  p  are  for  the  gaz  and  particulate  phase  in  the  two-phase  approximation, 
and  Cm  is  the  condensed  phase  /  gas  phase  mass  flow  rate  ratio. 


RECENT  DEVELOPMENTS:  FLUID-STRUCTURE  INTERACTION  [16] 

Fully  coupled  solution  of  fluid  flows  with  structural  interactions  is  called  fluid-structure  interaction. 
The  range  of  its  applications  is  important  in  many  engineering  disciplines  [17].  Some  current  applications 
are  pressure  waves  in  a  piping  system,  sound  waves  traveling  through  fluid-solid  media,  biomedical 
problems  such  as  blood  flow  in  a  diseased  artery  or  coupled  instabilities  in  power  systems.  Computer- 
aided  techniques  for  design  optimization  have  been  much  promoted  over  the  past  decades  and  have 
subsequently  reached  a  high  level  of  sophistication  within  many  single  disciplines  such  as  fluid  or 
structural  mechanics.  Flowever,  because  of  complexity  and  computational  cost  issues,  most  often  the 
coupling  effects  are  neglected.  For  example,  in  aerodynamics  optimization,  the  structure  is  assumed  to  be 
rigid.  Another  reason  for  separate  treatment  is  the  schismatic  split-up  of  engineering  disciplines, 
which  makes  it  difficult  for  one  person  to  have  in-depth  knowledge  on  all  of  these.  With  the  computer 
power  increasing  and  the  advent  of  parallel  processing,  research  on  fluid-structure  interaction  in  the  field 
of  numerical  aeroelastic  simulations  has  received  growing  interest  in  the  last  ten  years. 

In  fluid-structure  interactions,  the  combined  effects  of  inertial,  elastic  and  aerodynamic  forces  impact  the 
movement  of  the  fluid  and  solid  boundaries.  The  fluid  movement  exerts  aerodynamic  forces  on  the 
structure  that  reacts  and  in  turn  forces  the  flow  to  evolve  at  the  interface  with  an  interface  velocity. 
This  produces  the  coupling  effect  and  suitable  computational  strategies  need  to  be  developed.  The  problem 
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of  the  motion  of  the  fluid-structure  interface  that  occurs  in  coupled  aeroelastic  problems  is  generally 
addressed  by  solving  the  fluid  equations  on  moving  dynamic  meshes  with  an  Arbitrary  Lagrangian 
Eulerian  formulation. 

Fluid  structure  computation  has  applications  in  solid  propellant  grain  mechanical  design  and  unsteady 
simulations. 

Main  Features  of  the  Modeling 

Mathematical  models  of  coupled  problems  are  usually  coupled  partial  differential  equations  in  space  and 
time.  Often,  different  discretization  techniques  are  used  on  the  different  physical  components.  The  most 
difficult  part  of  handling  numerically  the  fluid-structure  coupling  arises  from  the  fact  that  the  structural 
equations  are  usually  formulated  with  lagrangian  coordinates  while  the  flow  equations  are  expressed  using 
eulerian  coordinates.  These  physically  heterogeneous  system  components  are  computationally  treated  as 
isolated  entities  that  are  separately  advanced  in  time.  Interaction  effects  are  viewed  as  forcing  effects  that 
are  communicated  between  the  individual  components. 

Arbitrary  Lagrangian  Eulerian  Formulation 

The  ALE  formulation  consists  in  solving  the  conservation  equation  on  a  moving  grid.  Since  if  the  grid  is 
fixed,  the  method  is  called  eulerian,  and  the  method  is  called  lagrangian  for  grid  points  having  the  material 
velocity,  the  ALE  formulation  is  a  generalization. 

Numerical  fluxes  must  be  computed  correctly  through  the  moving  faces  [17]. 

In  order  to  solve  the  problem,  the  computational  mesh  has  to  be  moved  or  deformed  during  the  time 
integration  of  the  fluid.  A  common  technique  to  deform  a  mesh  is  the  spring  analogy.  The  force  exerted  by 
the  nodes  j  which  are  connected  to  the  node  i  is  mathematically  expressed  by: 

j 

where  Ktj  is  the  stiffness  of  the  spring  between  nodes  i  and  j.  At  equilibrium  state,  the  force  at  every  node 
i  has  to  be  zero.  After  regrouping  the  terms,  the  iterative  equation  to  be  solved  yields: 

/2>. 

j  /  j 


Structural  Model 

The  structural  model  can  go  from  a  simple  reduction  of  stiffness,  mass  and  damping  matrices  of  the  finite 
element  mechanical  model  on  the  fluid  boundary  to  the  full  resolution  with  a  non  linear  structural 
mechanics  code. 

Fluid-Structure  Coupling  Algorithm 

The  coupling  algorithm  is  performed  in  a  staggered  way.  Structure  and  fluid  are  integrated  on  the  same 
time  scale  (given  by  the  fluid  time  step  required  by  the  fluid  solver)  although  the  characteristic  scales  for 
the  structural  system  and  for  the  fluid  flow  may  be  very  different. 
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The  fluid  and  structure  equations  are  coupled  by  imposing  that  at  the  boundary,  the  structure  stress  tensor 
is  in  equilibrium  with  the  fluid  pressure  and  that  wall  boundary  condition  occurs  on  the  fluid-structure 
interface.  It  implies  that  the  forces  and  energies  exchanged  at  the  fluid- structure  interface  are  balanced. 
Moreover,  mesh  and  structure  motions  are  also  coupled  by  continuity  conditions  on  the  interface. 

A  complete  cycle  takes  place  as  follows: 

1)  The  fluid  transmits  to  the  structure  a  pressure  profile  obtained  at  time  n  to  evaluate  the  pressure 
forces  exerted  by  the  fluid  F" . 

2)  The  structure  model  determines  the  displacements  q  "+1 . 

3)  This  structure  configuration  is  transmitted  to  the  fluid  model. 

4)  Fluid  variables  W  are  then  advanced  at  time  n  + 1 . 

5)  Back  to  step  1. 

The  position  of  the  structure  at  time  n  +  1  is  advanced  with  pressure  force  input  F"  defined  from  a  fluid 
pressure  profile  and  then  the  flow  state  vector  W"+l  is  computed  from  the  mesh  configurations  x"  =  q" 
to  xn+1  =  q""  .  The  position  of  the  dynamic  fluid  mesh  does  not  lag  behind  that  of  the  surface  of  the 
structure. 

With  no  damping  matrix  \D\,  the  structural  energy  expresses  as: 

E,=Fq{K\q  +  Fq{M\q+'qF, 

and  its  variation  during  a  time  step  is  AEs  =  E"+l  —  E"  =  {'  qn,i  —  q" )  F"  deriving  from  the  structural 
time  integrator.  On  the  other  hand,  the  transferred  energy  through  an  element  of  the  fluid-structure 
interface  can  be  written  as  AE,  =  ('x" 1 1  x" )  P  where  P  is  the  nodal  fluid  force  whose  expression 

depends  on  the  fluid  pressure  values  used  by  the  flow  solver  to  compute  the  fluxes  across  the  fluid- 
structure  interface.  The  force  and  energy  exchanged  between  the  fluid  and  the  structure  at  their  interface 
must  be  opposed.  To  ensure  these  principles,  a  procedure  enforcing  momentum  or  energy  conservation  is 
used  at  the  interface  in  step  4.  To  do  it,  fluid  pressure  involved  in  the  boundary  flux  is  re-adjusted. 
The  simplest  way  is  to  estimate  P  and  F"  with  the  mean  value  of  the  gas  pressure  available  at  the 
beginning  of  cycle. 

CONCLUDING  REMARKS 

In  this  lecture,  on  overview  of  general  problems  and  models  used  in  CFD  code  for  SRM  steady  interior 
flows  has  been  made.  The  second  paper  will  focus  on  unsteady  phenomena. 
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